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Abstract: An extreme mass ratio inspiral consists of two parts: adiabatic inspiral and 

plunge. The plunge trajectory from the innermost stable circular orbit (ISCO) is special 

(somewhat independent of initial conditions). We write an expression for its solution in 

closed-form and for the emitted waveform. In particular we extract an expression for the 

0^ associated black- hole ringdown amplitudes, and evaluate them numerically. 
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1. Introduction and summary 

Gravitational wave observatories (see for example the reviews [1] and references therein) 
demand knowledge of the waveform emitted by an inspiralling binary system of compact 
objects. Analytic control is possible in two limits. The first is the Post-Newtonian (PN) 
approximation, which holds whenever the velocities are small compared with the speed of 
light. The second limit, the one relevant to this paper, is that of an extreme mass ratio 
(EMR), for example, a compact projectile such as a stellar size black hole (BH) falling into a 
super-massive black hole. Numerical solutions complement the analytic methods and have 
become possible in recent years - see the review [2] and references therein. However, while 
numerical simulations work very well for comparable masses, they encounter problems in 
the EMR case due to the existence of two disparate scales which proves problematic for the 
discretized grid. So it happens that just where the numerical method loses precision, the 
analytic EMR approximation becomes more effective. Moreover, analytic methods improve 
insight into the problem, especially into its dependence on parameters. 

The evolution of an Extreme Mass Ratio Ispiral (EMRI) is customarily divided into 
two stages: an adiabatic inspiral where the system moves on quasi-bound orbits and slowly 
loses energy to gravitational waves, and a plunge phase where the system is set on a course 
of collision even when the self-force is neglected. In this limit the plunge phase simplifies to 
consist of geodesic motion of the projectile till it reaches the horizon of the larger BH. The 
ISCO (innermost stable circular orbit) is the border between the two stages, after which 
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the compact object freely falls into the horizon, followed by merger and ringdown (through 
quasi-normal modes) into a new stationary state. The smooth transition between adiabatic 
inspiral and plunge was beautifully understood in [3]. While the adiabatic inspiral is well 
within the applicability domain of the heavily studied PN approximation, known today up 
to order 3.5PN (as summarized in the review [4]) considerably less is known analytically 
about the plunge phase which is the subject of this paper. 

We would like to have a mathematical description of the waveform throughout its 
evolution, beginning with the adiabatic inspiral and lasting until the waveform vanishes 
after the black hole rings down and reaches steady state. 

The current state-of-the-art for such an all-time analytic expression goes under the 
name "Effective One Body formalism" (EOB) [5, 6], see also the recent review [7]. The 
EOB approach combines quite a number of ideas and ingredients, but here it will suffice 
to concentrate on its handling of the ringdown phase, which is quite independent of the 
other ingredients, and is based on the following assumption: the modeled waveforms during 
the adiabatic and the ringdown phases may be matched together without the need for an 
independent modeling of an intermediate phase. This assumption was later vindicated 
by full numerical simulations. More specifically, the adiabatic phase is carefully modeled 
and a general expression for the ringdown phase is matched directly onto it. In [8] the 
full waveform was computed in the extreme mass ratio limit by numerically evolving the 
trajectory within the EOB approach and numerically determining the radiation. In [9] the 
model evolved to its present form where the waveform emitted during ringdown is modeled 
by the most general ringdown signal - a linear superposition of its quasi-normal modes (eq. 
(18) there) 

^ gd °™W =J2 R n e " CT " ' + E R n ' . (1-1) 

n n 

where = a n ± ioj n are positive/negative quasi- normal-mode (QNM) complex frequen- 
cies, are the corresponding ringdown amplitudes, and n labels the QNM overtone 
number. This matching is implemented by matching the EOB model for the adiabatic 
phase with a general ringdown waveform (1.1) 2 around some judiciously chosen matching 
time, 3 matching at as many time points as the number of amplitudes to be determined. 

In this paper we shall compute the ringdown amplitudes directly from the theory for a 
certain special plunge trajectory. 4 For simplicity we consider the trajectory to be moving 
in the background of a non-rotating (Schwarzschild) BH, and we moreover consider the 

1 The method's main idea is an attempt to carry over to General Relativity (GR) the successful classical 
reduction of the two body problem to a reduced one body problem, where here in GR the reduced point 
particle moves in some "corrected" metric background. This background does not solve Einstein's equations 
and as such seems to offer enough adjustable parameters to fit for the sought waveforms. Its proven success 
makes it undisputable as an ad-hoc method. Yet it is not obvious that it is theoretically justified or natural, 
in the sense that a natural approximation would not contain unnecessary adjustable parameters nor would 
it lack parameters for higher order approximations. 

2 With R- = 0. 

3 The crossing of the "light ring" at r — 1.5r s . 

4 A similar approach was taken in [10]. However they differ significantly from the current study both in 
focus and in results. In particular, they do not calculate the ringdown amplitudes. 
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plunge trajectory to be the one which starts at the innermost stable circular orbit (ISCO) 
a t nsco = 3r s (in Schwarzschild coordinates, where r s is the Schwarzschild radius). This 
trajectory is special since the eccentricity of an orbit is known to decrease during the 
Keplerian regime r » r s of the inspiral [11], thus being a "Keplerian attractor" and 
hence only weakly dependent on initial conditions. In the case that the inspiral has been 
going on for enough time for the eccentricity to be essentially radiated away by the time 
ISCO is reached, then the plunge trajectory will indeed start from roughly ISCO. However, 
there are relatively recent indications that in some probable astrophysical scenarios other 
initial conditions are expected, namely initial conditions such that the last stable orbit is 
not circular, see for example [12, 13, 14]. Yet, in this paper we confine ourselves to the 
post-ISCO plunge for its simplicity and for being a Keplerian attractor. 

For the post-ISCO plunge trajectory we compute the form of the emitted gravitational 
waves. In the t — >■ oo limit we may extract the post-ISCO ringdown amplitudes, while 
the formal limit t — > — oo corresponds to radiation from ISCO. We wish to stress that 
the waveform under study continuously interpolates between two time regions: in early 
time the compact object is near ISCO and the radiation is characterized by the frequency 
of that trajectory, in the transition period the trajectory is not periodic and hence there 
are no sharp frequencies in the waveform, while at late times the decaying radiation is 
characterized by the (complex) frequencies of the quasi-normal modes which are set by the 
black hole, and not by the trajectory. 

In section 2 we solve analytically for the plunge trajectory, namely the trajectory with 
energy and angular momentum identical to the ISCO values. The expression for r = r((f>) 
is known (at least since [15]) to be especially simple and the time dependence can also 
be expressed in closed form (2.10,2.12). In section 3 we review the Zerilli/Regge- Wheeler 
[16, 17] theory for radiation in the background of a Schwarzschild BH for a given source 
of energy- momentum (using Martel-Poisson [18]). We proceed to substitute in the energy- 
momentum for our post-ISCO plunge trajectory. In section 4 we write a formal solution to 
the wave equation in terms of Green's functions (4.12). We extract the late-time ringdown 
form by deforming a frequency integral into the complex domain and transforming it to 
a sum over residues at the quasi-normal modes. The final expression for each post-ISCO 
ringdown amplitude (4.18) calls for solving certain ordinary differential equations (the 
radial wave equations) at the QNM frequencies and performing a certain weighted radial 
integral over the source. In section 5 we numerically evaluate these integrals and obtain 
the amplitudes of the leading ringdown modes shown in tables 1,2 and in figure 3. These 
values were confirmed by full numerical simulations [30]. 

Our expression (4.12) contains more than the ringdown amplitudes: it describes the 
full plunge waveform as a function of time. Solving for all (I, m) is equivalent to solving 
an inhomogeneous 3+1 wave equation and extracting the asymptotic outgoing signal. 

The impact of our work is that hereafter the EMRI post-ISCO ringdown amplitudes 
should be considered to be determined and known (at least under the assumptions above) 
and any waveform model must conform with it at t — > oo (and even as early as ISCO is 
crossed assuming the full plunge waveform is used rather than only its late time part). 

Generalizations. Several generalizations remain. First, one may consider the back- 
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ground to be a rotating BH (Kerr), in which case there are two kinds of post-ISCO plunge 
trajectories: co-rotating and counter-rotating, and the Zerilli/Regge- Wheeler equations 
must be replaced by the Teukolsky equation. Second, one may consider orbits where ec- 
centricity (and being off the equator plane in the Kerr background) was not washed away 
yet. Finally higher order effects may be incorporated such as accounting for the projectile's 
spin. 

2. Plunge trajectory 

Consider an inspiralling binary system in the Extreme Mass Ratio Inspiral (EMRI) limit, 
namely rh <C M where m is the mass of the compact object falling into a BH with mass 
M, or Schwarzschild radius r s = 2GM. 5 In order to calculate the gravitational radiation 
emitted from a system we need to find the dynamics of the compact source. For reasons to 
be discussed below we consider the (geodesic) trajectory which spirals out of the innermost 
stable circular orbit (ISCO) at nsco = 3r s and plunges into the black hole. In this section 
we shall solve for this trajectory. 

We concentrate on the case when the inspiral has been going on long enough for 
the radiation reaction force to have circularized the trajectory by the time the compact 
object reaches the innermost stable circular orbit (ISCO) at risco = 3r s . This assumption 
deserves a discussion. Peters [11] found that in the Keplerian regime r 3> nsco energy 
and angular momentum are lost such that eccentricity decreases with the following rate 

12 / 121 n\ 2299 

p = c eT9 ^l + _ e 2 J (2.1) 

where p is the Keplerian semi-latus rectum, and e is the eccentricity, such that the trajectory 
is given by r{9) = p/ (1 + ecos(#)) and the semi-major axis a is given by a = p/(l — e 2 ). 
Assuming small e and extrapolating (2.1) down to ISCO we can estimate that 



19 



e/ ~ ei 



(rfc>) » ^ (22) 

where /, i subscripts denote final and initial, respectively. Note that this is only an approxi- 
mation since the radiation reaction force is different in the strong gravity regime r > nsco > 
and it is actually known that e increases somewhat just before ISCO [20]. Altogether, (2.1) 
and (2.2) quantify the rate of circularization and explain why e would be small at ISCO 
given pi 3> nsco- Hence ISCO is a Keplerian attractor. Yet, as mentioned in the intro- 
duction it should be borne in mind that initial conditions such that the last stable orbit is 
not circular are also to be expected. 

Since ISCO is a (marginally) unstable orbit, we choose the plunging orbit with the 
same conserved parameters as those at ISCO, namely 

E = Ejsco = 

L = Lisco = V3r s , (2.3) 



J We reserve the letter m to denote the magnetic number - an index for the spherical harmonic functions 

-ylm 
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where -Eisco> -^isco are the particle's energy per unit mass and angular momentum per 
unit mass at ISCO. During the plunge radiation reaction can be neglected as the object is 
already falling "without its help". 6 

We wish to find explicitly the trajectory of this free-fall. Without loss of generality 
the trajectory can be considered to lie in the equatorial plane and the equations governing 
the object's motion are 

E = f(r) £ (2.4) 

(2.5) 

E 2 = (^) 2 + f(r)(P/r 2 + l), (2.6) 

where f(r) = 1 — ^f. 

In our special orbit, there is a simple relation between the r and <j) coordinates of the 
particle (noticed at least as early as [15]) 

™ = 1 + <*^F- (2 - 7 > 

where <j>Q is an arbitrary constant (sometimes we will choose it to be zero). We will also 
use conventions such that ■£ > 0. 

We proceed to solve for the relation between t and r through integration of equation 
(2.6), using equation (2.4) to change the derivative with respect to r to a derivative with 
respect to t. For a generic orbit we can solve for t = t(r) by quadrature involving an elliptic 
integral 

E I dr = - f dt . (21 

/(r)^p-/(r)(g + l^ 

For our special orbit, the expressions simplify due to the special form of the effective 
potential (namely that both its first and second derivatives vanish at nsco)- Plugging in 
our constants of motion, -Eisco an d -^iscO) it can be seen that (2.6) becomes 



This obviously simplifies the integrand, and with a change of variables 
X := 2 ( ris f° — l) it is brought to the form of a rational function. Evaluating the integral, 

6 The leading corrections for the ISCO values at plunge were computed in [3] eq. (3.26) and they were 
found to be proportional to (m/M) 4//j . 
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we find that our special orbit satisfies the following relations 



t(r)/r s 



2r ( 1 — 12 ^0 / \ 

i r ' - 22V2 tan" 1 ( y^) + 2 tanh" 1 (y/x) + t (2.10) 



r s JX 



mir. = 3V6 (0 - M (l + (<t _^, + 12 ) - ^ tan- (^\ (2.11) 



+ 2 tanlT 1 f , — , ) + i 



(2.12) 



where to is an arbitrary constant which we shall fix in section 5. We note that x ls closely 
related to the standard change of variables in Kepler's problem u := 1/r and that x varies 
from x = at ISCO to x = 1 a ^ the horizon. 



3. Source for gravitational waves 

In this section we proceed to describe the equation satisfied by the emitted gravitational 
waves given a source moving on the geodesic trajectory described in the previous section. 

3.1 A general source in the Schwarzschild background 

The theory of linearized perturbations of Schwarzschild space-time was initially developed 
in the pioneering works of Regge and Wheeler (1957) [16] and of Zerilli (1970) [17]. Essen- 
tially, when expanding the metric perturbations in spherical harmonics, the physics of the 



(a) (b) 




Figure 1: (a) Spatial trajectory of the plunging particle, r((f>), normalized so that r s = 1. The 
dashed line at r = 3 represents the ISCO, whereas the dotted line at r = 1 represents the horizon, 
(b) r coordinate of the plunging particle as a function of coordinate time t (equation (2.10)). 
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perturbations reduces to the physics of a (1 + 1) dimensional wave equation with a nontriv- 
ial potential, satisfied by two "master functions" (for odd and even parities) constructed 
from the metric. These master functions characterize the metric perturbations. 

A convenient formalism, suitable for the problem of analyzing gravitational waves gen- 
erated by material sources, was introduced by Martel and Poisson (2005) [18] . That paper 
presents gauge-invariant and covariant form of the Regge- Wheeler equation (describing 
odd-parity perturbations of Schwarzschild space-time) and of the Zerilli equation (describ- 
ing the even parity perturbations) with two corresponding scalar master functions and 
invariant source terms constructed from the stress-energy tensor of the matter responsible 
for the perturbation. In this subsection we review the general theory of Schwarzschild BH 
perturbations with sources, closely following [18] (see also [19]). 

The Schwarzschild metric can be written as 

ds 2 = g^ dx^dx u = g ab dx a dx b + r 2 Q AB d6 A d9 B (3.1) 
g ab dx a dx b = -f dt 2 + f' 1 dr 2 (3.2) 
fl A B d9 A d6 B = d9 2 + sin 2 9 d<p 2 , (3.3) 

where upper-case latin indices run over the values 2,3 (9,<j>) and lower-case latin indices 
run over the values 0,1 (t,r). We study small perturbations around the Schwarzschild 
space-time, so it is useful to write the metric in the form g^ v = g^ v + h^ u , {h^l <C 1. Next, 
we expand the perturbation in spherical harmonics, according to the transformation laws 
of the different components under rotation. For the odd-parity sector 

Kb = (3.4) 
h*B = Y,ti? xl B ( 3 - 5 ) 

Im 

h AB = J2h l 2 m X%, (3.6) 

ill) 

where X 1 ^ 1 and X 1 ^ are the odd-parity vector and tensor spherical harmonics 7 , respec- 
tively. For the even-parity sector 

fc«6 = 2>S? ylm (3.7) 

Im 

KB=J2^ Y B m ( 3 - 8 ) 
Im 

h AB =r 2 Y, ( Rlm ^AB Y lm + G lm Y l A % , (3.9) 

Im 

where Y lm , Y B m and Yjfg are the even-parity scalar, vector and tensor spherical harmonics, 
respectively. Plugging these expressions into the Einstein equations and using the orthog- 
onality relations of the spherical harmonics, a wave equation for every (l,m) is obtained, 



7 More details regarding the spherical harmonics can be found in [18], and in particular we record their 
normalization conventions: §Yi m Y Vm > dSl = 5 u >5 mm t, j Y^YX" 1 dD. = 1(1 + l)5 u >5 mm t, j X^X^" 1 dtt = 
1(1 + 1) 5 u ,5 mm ,, fYj^YXF <Kl = \(l-l)l(l + l)(l + 2)5 nl 5 mm , and finally / X^X 1 ^' dO. = |(Z- 1)1(1 + 

' u mm' • 
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satisfied by the corresponding master functions. The Cunningham-Price-Moncrief function, 
characterizing odd-parity perturbations, is defined by 



lm ._ 



1r / 2 

- ab (d a h l b m -- 
\ r 



(3.10) 



(I -!)(* + 2) 

where e ab is the Levi-Civita tensor on the (t, r) manifold and r a = (0, 1) is the unit co- vector 
in the r direction. The Cunningham-Price-Moncrief function satisfies the Regge- Wheeler 
equation 

(n-vt)ti m = s l m , (3.ii) 

where 



1(1 + 1) 3r s 



(3.12) 



is the Regge- Wheeler potential, and S l ™ is the source term, on which we will focus later. 
The Zerilli-Moncrief function, characterizing even-parity perturbations, is defined by 



lm 



2r 



1(1 + 1) . 



K lm + j(r a r b (h% 



2V„ - r 2 V,,C 



r r 



'V a K lm ) 



where 



A:= (l-l)(l + 2) + 



3r, 



(3.13) 
(3.14) 

2 -v- ■ r „■ , , (3-15) 

and r a = (0, 1) is the unit co-vector in the r direction on the (t, r) submanifold. The 
Zerilli-Moncrief function satisfies the Zerilli equation 



A ":- K+\l(l + l)G-\r a [j a -^r 2 \ ll G 



where 



Vo 



1 

A2 



At + 2 3rs 



+ 



9H 



(3.16) 
(3.17) 



is the Zerilli potential (/j, := (I — l)(l + 2)) and S^ 7 ™ is the source term. 

The source terms, as was earlier stated, are constructed from the stress-energy tensor 
of the material source of gravitational waves. In the odd-parity sector, the source term is 
given by 



lm 



where 



S, 



pa 



2r 



(*-!)(* + 2) 

167TT 2 



e ab V a P b , 



1(1 + 



-I 



T aB x *lm dQ 



(3.18) 
(3.19) 



and X* denotes the complex conjugate of X. 

In the even-parity sector, the source term is given by 



12r, 



-2 r 2 r a V a Q + r a r b Q ab + 2 r / Q b 



A 



+ 



At (ax - 2) + 6 (fi - 3) j + 21 ^| 



21 



(3.20) 
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where 



Q ab = 8vr j T ab Y* lm dn (3.21) 

Q a = 1 ^ 7rr ' 2 j T aB Yg m dU (3.22) 

Q b = 8vr r 2 J T AB Q AB Y* lm dQ (3.23) 

^ = {l -rZm + 2) f TABY ^ in (3 ' 24) 

Q = 9 ab Qab ■ (3.25) 



In this work we ran a check, independent of the derivation of results in [18], of the 
expressions for the source terms (3.18) and (3.20). The idea is simple - we checked that 
both sides of the Regge- Wheeler and Zerilli equations with source terms, equations (3.11) 
and (3.16), agree - when written explicitly as functions of the coefficients of the spherical 
harmonics in the decomposition of the metric perturbations (e.g. h 1 ™,...) appearing in 
equations (3.4) - (3.9). When substituting expressions (3.10) and (3.13) in the left hand 
sides of the equations (respectively) such an explicit expression is obtained. For the right 
hand sides of (3.11), (3.16), we used relations between the coefficients in the spherical 
harmonic decomposition of the stress-energy tensor (P a in (3.19), the Q's in (3.21)-(3.24)) 
and the metric perturbation coefficients (h l ™, etc.). These relations can be obtained 
by substituting the decompositions to spherical harmonics in the Einstein equation and 
equating the different spherical harmonic coefficients. Although the idea is simple, its 
execution demands some tedious algebra. The process described above was carried out for 
equations (3.11), (3.16), and the two sides of both equations were found to coincide. 

Combining the two wave equations (3.11), (3.16) into a single notation and making 
the wave operator explicit, we obtain in the frequency domain 

V#.(r,u) = /^7.(r,a;), (3-26) 

where r* = r + r s log(r/r s — 1) is the usual tortoise coordinate and ^^(r, cj) axid S l J?(r,cj) 
are the Fourier transforms of the master function and source term, respectively. Notice 
that changing to derivative with respect to r* instead of r casts the equation explicitly in 
the form of a flat space 1 + 1 dimensional wave equation with a potential. The integration 
of these equations will be discussed in section 4. 

For later use we record the relation between the master functions and the transverse- 
traceless part of the metric perturbation (which is gauge-invariant) at future null infinity 



Kb = (3.27) 

h aB = (3.28) 

h A B = rJ2 (^Yab + ^X%) . (3.29) 

Im 
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3.2 Plunging mass as a source 

We turn to calculate the source terms corresponding to our plunging point particle. The 
stress-energy tensor corresponding to a point-like mass is given by 

T^ = rh J dr (-g)~ l/2 u^u" <5 (4) (x a - x°(t)) , (3.30) 

where 5^(x a — y a ) is the four-dimensional Dirac functional, is the 4- velocity and x p (r) 
is the location of the point particle which approximates the compact object. x p (r) are the 
same functions discussed in section 2 and the added subscript p is to distinguish them in 
the current context. 

In the Schwarzschild space-time, after localizing the r integral onto t p (r) 

T^ = rhL U " {t) f {t) 5 (r - r p (t)) S {d - |) 5 (</> - fr{t)) . (3.31) 

To compute the odd-parity source term, we first need to calculate P a (3.19). The integral 
over the angles is trivial due to the delta functions and we obtain 

pa = l ^™L f{r P ) ^ s (p _ rp{t)) m m) e - imMt) (3 32) 

l[L + L) hi r~ 

where B{l,m) is defined by 

= B{l,m)e- im,j> , (3.33) 



\r*lm 



and can be explicitly written as 



s 21 + 1(1 -m)\ (l + m)\\ /7T,, 



for Z > m (for Z = m replace the factor of (I — m — 1)!! in the denominator with \/tt). 

Equation (3.18) gives us an expression for the odd-parity source term as a function of 
r and t 

5 °><*) = (j _ !)( f + 2) [r^t^ + ^C/P*)] • (3.35) 
The source term in the frequency domain is given by 

ylm f„ , ,\ / iwt olm 



/oo 
e-™Sfy e (r,t)dt . (3.36) 
-oo 



To deal with the time derivative in the integrand of (3.36), coming from the first term in 
(3.35), we integrate by parts to get a factor of — iui from the exponent's time derivative. 
Now we perform a change of variables in the integral, t — > r p (t). Additionally, we have 
explicit analytical expressions for the functions t p = t p (r) (equation (2.10)), appearing in 
the integrand in the factor e -* aJ *p( r ) ) an d ^ p = <j) p {r) (equation (2.7)), appearing in the 
factor e - 4m< M*). Using these relations, all the dependence of the integrand on r p (t) is 
explicit. Now, using the delta functions, the integral can be evaluated analytically. 
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The final result for the frequency domain odd-parity source term is, for r s < r < 3r 
32^/3 7rmr, 



Ar 



B(l, m) e 



-i(ui t p {r)+m <j>p{r)) 



-IU) 



r 

3r s — r 



1 + 3 



(3.37) 

6\/6 r s r y/2 r 1 / 2 (3r s - 4r) 
(3r s — r) 3 (3r s — r) 5 / 2 



where A = (Z - 1)1(1 + 1)(Z + 2). For r g (r s ,3r s ), S l m (r,co) = 0. 

The calculation of S'g m (r, w) is done in a very similar way (but it is a bit longer). 
Instead of calculating P a we will need to calculate the Q's from equations (3.21)-(3.24). 
We will also write explicit expressions for the relevant spherical harmonics in the 9 = | 
plane 



Y 



*lm 



Y 



Y. 



*lm 



4m 



A(l,m)e- im * 
-imA(l,m) e- im4 > 



-m 2 + -1(1 + 1)) A(l,m)e- i 



where 



A(l,m) := 



'2 +1 (l-m ) (Z + m-1) /7r„ 

— r^T^ fr St TTi — sin ( -( + m + 1 

4vr (l + m)\ (l-m)\\ V2 V ; 



(3.38) 
(3.39) 

(3.40) 
(3.41) 



After calculating the time domain source term by equation (3.20), we Fourier transform 
it, in the same spirit as in the odd-parity case (using the same change of variables), to obtain 
an explicit analytical expression. 

The final result is, for r s < r < 3r s 



S l e m (r,u>) 



167r rh 



^A(l,m) e 



-i{u!t p {r)+m<f> p (r)) 



(ji + 2)Ar 

9rr 2 A (-2m 2 + ll + 2) 6i^/¥ (r 2 + 3r 2 ) (3y/3mr a (r, - r) - 2 v / 2r 3 «;) 



ii(3r s — r) 3 / 2 



(3r s — r) 3 



WW*. + 3 ( '' 2 + 3r ' 2) 2) " + 6 :i- ( " - 3) + 21ra (3.42) 

rA(3r s — r) 3 / 2 

3(r - r s ) (4r 3 - 3r 2 r s + 24rr 2 - 45r 3 ) + 18r 2 (r s - r) 4r s (3r s - r) 3 / 2 



(3r s — r) 5 / 2 



(3r s — r) 3 / 2 



rA 



For r £ (r s ,3r s ), S l c m (r,uj) = 0. 

Note that S e is non-zero only for even I + m and similarly for S . 



4. Solving the wave equation 



In this section we discuss the solution to the wave equation (3.26). 



4.1 Green's function 



Knowing the explicit form of the source terms we proceed to calculate the gravitational 
waveform. The boundary conditions for equation (3.26) are outgoing waves at infinity 



and ingoing boundary conditions at the horizon 

The solution to equation (3.26) is given by 

^ e (r,u>) = [ G l 0/c {u,r^)fS l J c {<Mdr: , 



oo 



-oo . 



(4.1) 



(4.2) 



(4.3) 



where G l Q , (u;,r*,r*) is the frequency domain Green's function given the above-mentioned 
boundary conditions. It satisfies 



Or. 



+ u 2 -fVl /c (r) 



G l o/c (uj,n,ri) = 5{n-ri) 



(4.4) 



with the same boundary conditions. 

In order to construct the Green's function we must first define two independent solu- 
tions to the homogeneous equation (c.f. [21]) 



dr 2 



+ u 2 -fVl /c (r) 



u l>/e( w » r *) = • 



(4.5) 



Let us denote by u 1 ^ the solution to (4.5) which satisfies (we will sometimes suppress the 
parity) 



u^(r*,o;) ->■ e lujr * ; r* ->■ +oo . 
At the horizon Uoq can be expanded as follows 

uL(r tlW )^< t ( W )e~-«+4(")e' 



i / * 



-oo , 



(4.6) 



(4.7) 



where A- m and A out are some u dependent complex coefficients. In a similar manner, we 
denote by u l hor the solution to (4.5) which satisfies 



«hor( r *> w ) 



-oo 



At infinity it can be expanded as 



; r* ->■ +oo . 



(4.8) 



(4.9) 



The Green's function satisfies (4.4) with outgoing wave boundary conditions at infinity. 
Therefore, for r* > the solution is proportional to u^. Similarly, for r* < the solution 
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is proportional to ui . What remains is to properly match the solutions at r* = r'^ and 
obtain the Green's function 

w ~ 1 (^) u Lv( r ^ U) ) ul oo( r '*,^) \r,<ri 

(4.10) 



where W(u) := u^dr^u 1 ^ — u^d^u^ is the Wronskian of u l hor , u 1 ^ and is independent of 
r*. Evaluating it both at the horizon and asymptotically we obtain VF(w) = —2icuBi n = 
—2ibjA out and in particular B{ n = A out . 

We are interested in calculating the gravitational waveform measured by a distant 
observer (r 3> r s ). We can, therefore, simplify the Green's function (4.10) 

G[ /e (u,,r.y.)~-—}— u{ or {r^) . (4.11) 

' 2,1 LO B in 

Knowing the Green's function (4.10), (4.11) we can proceed to express the solution to 
the wave equation. Using (4.3) and transforming back to the time domain we have 



1 



<TeM) = 2^/ dr> I e^G l o/c (co,n,ri)S l l(r:,co)du 



3/v 



(4.12) 



^iuit /~<l i , . „ J \ aim 

-oo 

4.2 Ringdown amplitudes 

It is known that as a perturbed black hole settles back to its stationary state the background 
metric exhibits an exponentially decaying radiation known as ringdown, or quasi-normal 
modes. This can be seen from the analytic structure of Green's functions for the Regge- 
Wheeler and Zerilli equations in the complex frequency domain. It is known (e.g. [22]) that 
these Green's functions have poles in the upper half plane (Im(w) > 0) and a branch cut 
along the positive part of the Im(w) axis as shown schematically in figure 2. The poles are 
interpreted as the quasi-normal modes. Their locations are the quasi-normal frequencies, 
whose imaginary part is responsible for their decaying nature at late times. The branch 
cut is interpreted as the "tail" resulting from scattering the waves off the background 
geometry at large radii. Here we shall proceed to derive a formula for the ringdown or 
quasi-normal mode amplitudes. Concentrating on them is justified since they are known 
to dominate over the tail for late times but not too late times - time region (iv) of [22], 
the discussion after eq. (46). See [23, 24, 25, 26] for further discussion. In addition, these 
amplitudes are one of the ingredients which go into current waveform models, as reviewed 
in the introduction, but as of now they are matched (with a different part of the modeled 
waveform) rather than computed from the theory. 

We now turn to evaluate the integral in square brackets in (4.12), utilizing the contin- 
uation to complex frequencies. We deform the contour from the real axis to parts (b) and 
(c) in figure 2 and collect the residue contribution from the poles. 

The poles correspond to zeros of the Wronskian. We denote by uj n i the ro-th pole for a 
given I (ordered by increasing imaginary part). At these frequencies B\ n {uj n i) = A out (uj n i) = 
and the leading behavior for B^ n is 

-Bin ^ du>Bm\ U =u M {U-Unl) ■ (4.13) 
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Im(fti) 



-.. (b) 



(c) 



Re((y) 



Figure 2: Analytic structure of the Green's function in the upper half ui plane and the integration 
contour. The thick line represents the branch cut on the positive Im(cj) axis. The dashed line 
represents the integration contour. The dots represent some of the Green's function's poles at the 
QNM frequencies. 



For convenience we assign notation to these constants 

Pnl '■= —dujBi, 



(4.14) 



Their formalism and numerical computation can be found at [27]. 

The r* 3> r s Green's function (4.11) (implying r* > r*) for oj ~ uj n i can now be written 

as 



G o/e (u,r*,r{) 



(4.15) 



2iu n i j3 n i(u} - UJ n i) 

At late times when the contribution of parts (b) vanishes and part (c) is neglected as 
we discussed above the integral over frequencies (4.12) can be approximated by the sum 
over the residues 

r-3r a 



1 

2tt 



dr'^TrzResU^ e*"* G l o/c (u>, U , S l %(C u) . (4.16) 



Plugging in (4.15), we obtain 



iu> ni (t-r*) /-3r 



2 /3 n / 



(4.17) 
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where the index n runs over all QNM's for a given I. 

We identify the coefficient of e tuJ m{t-r*) w ith the ringdown amplitude 



1 



3 7. 5 



Rnlm := ~ 77" / <oM) S$(K, Unl ) dr' . (4.18) 

* Unl Pnl Jr s 

This is our main result for the ringdown amplitudes. 

For r — > risco ^hor is finite and the third order poles in the source S are suppressed 
by the term exp(— iut) which is bounded by | exp(— iu)t)\ < exp(Im(w)t) ~ exp[Im(cj) • 
const/ \/r — nscoj an d so the integrand vanishes exponentially in the limit. 

In the limit r — > rh or , however, the wavefunction diverges 



f \ r s 



Uhor expfwjr*) oc 1 . (4.19) 

Vs J 

The source term, too, diverges in this limit in the same fashion. At r — > r s the trajectory 
satisfies r* ~ t + const and so 

1 I . (4.20) 

r -i \ r s 



Altogether, then, the integrand is oc — 1J s and so the integral diverges when Im(^-) > 

2- This spurious divergence needs to be regularized. This problem of a diverging QNM 
excitation integral was realized already in [22] for different excitation scenarios. Two differ- 
ent methods (equivalent in principle) to regularize the integral are discussed there, one of 
them being the divergence subtraction method (introduced in [29] for excitation integrals 
of non-QNM frequencies). This is the method we used in practice when evaluating the 
amplitudes numerically, and it will be discussed in the next section. 



5. Numerical Evaluation 

In this section we numerically evaluate the amplitudes (4.18). 
5.1 Method 

We started by calculating the numerical values of the QNM frequencies co n [. One can use 
either one of several methods (see the review [28]), and we chose Leaver's continued fraction 
method (see [22]) which can be extended to yield also the QNM wavefunctions Uh OT (u) n [, r). 
More specifically, the wavefunctions are obtained in the regime of interest as a power- 
series expansion in (1 — ) multiplying a factored-out non-analytical part. The calculated 
frequencies were confirmed by the literature, and the wavefunctions were confirmed by 
substitution into the homogeneous Regge- Wheeler and Zerilli equations. We continued to 
calculate the excitation factors (or transmission residues) f3 n i, again using Leaver's method, 
and the first few n's and Vs were confirmed by Leaver's results. We found [27] to be a very 
useful and detailed guide to the calculation of the excitation factors (there the more general 
rotating BH case is considered). 
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The next ingredient in the integrand is the source term S l ™ e (r*,uj n i), which is given 
analytically by (3.38), (3.43) where one should substitute the trajectory functions t p (r) 
(2.10) and <f> p (r) (2.7). 

It remains to perform the numerical integration in (4.18). As noted in the end of last 
section, this integral diverges and must be regularized by what is essentially an analytic 
continuation. To do so, we used the so-called divergence subtraction method ([29], [22]). 
The idea is the following: we set the integrand to zero at the horizon (r — > rh or ) by 
subtracting from it a function / that behaves the same as the integrand as r — > 
thereby getting rid of the unphysical divergence at that integration limit. In order not to 
add anything finite to the integral we will need to add the value of F := j fdr at the upper 
limit of integration. The integration constant in F is determined by the condition that F 
does not have a part ~ const, as r — > 1. That is, if / is of the form given below in (5.2), we 
integrate it using the rule J r a dr = without the addition of a constant. In summary 



3;- f 



physical part 



3n 



dr' + F\ 3rs 



(5.1) 

/ is, of course, not unique. Any function which will be easy to integrate analytically 
(and behaves in the desirable way as r — > r^or) will do. The (natural) / we chose to use 
was 

f=(--l) S (A + (A 1 (r-r s ) + ...) (5.2) 

where the coefficients A{ are determined from our expressions for it^ and S^? e and enough 
terms must be taken as to assure convergence. 



log 10 |/?llm| 
0.0^^ 



Figure 3: Ringdown amplitudes from tables 1,2 depicted on a logarithmic scale for leading n = 1, 
2 < I < 6 and several m values. The amplitudes decrease as I — m increases, for instance for I = 3 
the top point corresponds to m ~ 3 while the bottom one to m = 0. 
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5.2 Results 

From equation (4.17) it is clear that the ringdown amplitudes are defined only up to a shift 
in the time coordinate; if we take t —> t + t' the amplitudes will transform as 

Rnlm — t Rnlm& lU)nlt ■ (5-3) 

We fix to, the integration constant in equation (2.10), somewhat arbitrarily such that 
r(t = 0) = 1.1 r s , which we observed to be shortly after the onset of ringdown [30], namely 
the region which is well described by (4.16). 

It is worthwhile to note that the decay constants (imaginary parts of the frequencies 
u n [) do not vary with m and vary only little with I while they vary strongly with overtone 
number n and hence (5.3) shifting to will not affect the ratio of magnitudes of amplitudes 
which differ only by their m, it will weakly affect different Z's and finally strongly affect 
the ratio for different n's. Note also that the real signal is gotten by adding the complex 
conjugate to sum of ringdown exponentials. This is equivalent to the doubling of QNM 
frequencies R(uj) — > —$t(uj) as in figure 2. 

Results for amplitudes of the first few dominant modes are displayed in table (1). In 
table (2) the contribution of higher Z's is displayed, up to I = 10. In table (3) numerical 
values of QNM frequencies are displayed for the first few (n,l). Figure 3 gives a graphic 
illustration of several amplitudes on a logarithmic scale. 

In figure 4 we display the ringdown waveform (calculated from the amplitudes through 
equations (4.16), (3.29)), as viewed by observers in different angular positions. Notice 

K 

0.5 r 




Figure 4: Ringdown waveforms summing the dominant (I, m, n) contributions. The (distant) 
observer is situated in the equatorial plane (9 = J), at: (a) <f> = 0, (b) <f> = \ , (c) <f> = tt, (d) 
<j> = The horizontal axis represents t — r in units of r s , while the vertical axis represents h + , 
the + polarized part of the gravitational wave normalized by r/rh. Notice the asymmetric nature 
of the ringdown waveform. 
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the highly asymmetric nature of the ringdown waveform - observers in different azimuthal 
positions see very different waveforms. This is due to the fact that the trajectory breaks 
azimuthal symmetry, by infalling at some specific angular position. Of course, waveforms 
for observers in any angular position can be easily calculated from the amplitudes, via 
equation (3.29). 
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Table 1: Numerical values of the ringdown amplitudes R n i m f° r I — 2, 3, 4, 5 , several values of the 
overtone number n and all values of m. Conventions: to (2-10) was taken such that the plunging 
object is at 1.1 r s at t = 0, the overall normalization is defined by (3.10,3.13), and the normalization 
of spherical harmonics is given in footnote 7. 

/ = 2 



m 


n = 1 


n = 2 


n = 3 


2 


-0.0985724 - 0.747787* 


-0.229354 + 0.428849* 


0.167484- 0.251937*' 


1 


-0.0210521 + 0.399297*' 


0.441304 - 0.31877i 


-0.489451 + 0.053276*' 





-0.0887841 + 0.0979244* 


0.303357 + 0.0416042*' 


-0.249012 - 0.314569*' 


-1 


0.0274099 + 0.00889306* 


-0.0324417 - 0.0903849*' 


-0.0867015 + 0.154488* 


-2 


(-0.735088 + 3.59504i) x 10~ 3 


0.0135078 - 0.0089118i 


-0.0394189 - 0.010878* 



/ = 3 



m 


n = 1 


n = 2 


n = 3 


3 


-0.0375476 + 0.141024* 


0.117107- 0.057957*' 


-0.0862334 + 0.00839625* 


2 


0.0336116 - 0.0585892*' 


-0.104903 + 0.00738472* 


0.0855561 + 0.0540601*' 


1 


0.0186212 - 0.0145793*' 


-0.0499562 - 0.0209467* 


0.0245107 + 0.0656791* 





(-7.36412 + 0.205584i) x 10~ 3 


0.0119392 + 0.0185584*' 


0.0108096 - 0.0349259*' 


-1 


(-1.32958 - 1.74202*') x 10~ 3 


-2.4849 + 7.94279*') x 10~ 3 


0.0160937 - 0.00679286* 


-2 


(-3.59688 + 3.99339*') x 10~ 4 


(2.43071 + 0.097078*') x 10~ 3 


(-3.98823 - 4.75782*') x 10~ 


-3 


(-5.2293 - 6.09408*0 x 10~ 5 


(-0.284844 + 4.19003^) x 10" 4 


(9.4908 - 8.54244*') x 10~ 4 



1 = 4 



m 



n 



1 



4 0.0318095 - 0.0308272*' 

3 -0.0165716 + 0.00919402*' 

2 (-6.86211 + 1.83095s) x 10~ 3 

1 (2.32275 + 0.593617*0 x 10~ 3 

(4.47817 + 5.89624*0 x 10" 4 

-1 (0.514994 - 2.12774*0 x 10~ 4 

-2 (6.24213 - 1.42382*) x 10" 5 

-3 (-0.482705 - 1.61408*') x 10~ 5 

-4 (2.7626 - 0.0853029*) x lO" 6 



n = l 

-0.0489266 - 0.00923324*' 
0.0289206 + 0.0167393*' 
0.0114492 + 0.0128208*' 
(-1.76906-7.17998*) x 10~ 3 
(1.10773 - 2.61122*') x 10~ 3 
(-9.22639 + 3.19073i) x 10~ 4 
(-2.49301 - 2.12194*') x 10~ 4 
(-5.20619 + 8.05728i) x 10~ 5 
(-1.26026 - 1.17755i) x !0 -5 



m n = 1 

~5 -0.0166349 + 0.00302628* 

4 (6.80123 + 0.509818i) x 10~ 3 

3 (2.61522 + 0.629434*0 x 10~ 3 

2 (-7.32533 - 5.82304*') x 10~ 4 

1 (-1.14018 - 2.79764*') x 10~ 4 

(-3.06403 + 8.96566*') x !0~ 5 
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Table 2: Numerical values of the ringdown amplitudes R n i m for I = 6 to 10, overtone number 
n = 1 and dominant m's. 



1 


771 = / 


m = I — 1 


6 


(6.73817 + 3.1156*) x 10~ 3 


(-2.21976- 1.78938*) x 10" 3 


7 


(-1.75245-3.12505*) x 10~ 3 


(0.356901 + 1.26881*) x 10" 3 


8 


(-0.229465 + 1.83699*) x 10~ 3 


(2.27002 - 6.14468*) x 10~ 4 


9 


(6.97897 - 7.26025*-) x 10~ 4 


(-2.87291 + 1.88928*0 x 10~ 4 


10 


(-4.03242 - 3.58438*') x 10~ 4 


(1.04411 + 1.44422*) x 10~ 4 



Table 3: Numerical values of some dominant QNM frequencies in units of r 



1 n = 1 n = 2 n = 3 

2 -0.74734 + 0.17792*- -0.69342 + 0.54783* -0.60210 + 0.95655*' 

3 -1.19889 + 0.18540* -1.16529 + 0.56259* -1.10337 + 0.95818*' 

4 -1.61836 + 0.18832*- -1.59326 + 0.56866* -1.54542 + 0.95981* 

5 -2.02459 + 0.18974*- 

6 -2.42402 + 0.19053* 
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